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ABSTRACT 


In this thesis we deal with the generation of bivariate 
random variables with identical marginal distributions and 
specified correlation coefficients. Many schemes have been 
put forward for different cases; they almost always involve 
computing the inverse probability distribution of the 
marginal distributions or exploiting special properties of 
the random variables in question. 

A simple mixture-truncation method is put forward here 
to generate bivariate distributions with any marginal 


distribution from univariate generators. 
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I. INTRODUCTION 


In many simulation applications, it 1s required to 
generate bivariate random variables which have identical 
marginal distribution. 

For example, in reliability problems, the assumption that 
two components have independent exponential failure times is 
very often unrealistic, and much effort has gone into deriving 
bivariate exponential random variables to handle these situa- 
tions [Gaver (1972), Olkin & Marshall (1967), Downton (1970)]. 
Now except in very specific physical situations it may be 
difficult to specify the complete bivariate distribution of 
life time of each component. However, it may be realistic 
to specify the marginal distributions and some measure of 
dependence (usually the correlation coefficient) between the 
life time of each component. In this kind of situation, we 
can use bivariate random vectors having given marginal dis- 
tribution and dependence to solve the problem in simulation. 
To generate these vectors, there exists some previous works 
and existing methods, discussed in Section II. As seen in 
Section II, most previous work is specific to specified 
marginal distributions and uses inverse transformation methods 
as a baSic concept. 

An example is the recent work by Johsnon and Tenenbein 
(1979), to generate a bivariate random vector (X,Y) which has 
marginal distribution Fi(x), F, Cy) and correlation p, by a 


weighted linear combination method. 





Define 


X = F)"(H, (U)) 
Y= F>*(H,(¥)) 

Or 
v= Pot (1 SEG 


where Hy and H. are Ener ocumulative distribution functions 


Wera.c.) Of U and V respectively and 


Wee —eehe + ( | ely 


where U', V' are i.i.d. random variables with probability 


=a -l] 
1? F, , and ¢ 


are specific to the marginal distribution and correlation 


Sem@erty functino g{( ). In this procedure, F 


desired. The functions a and FS are difficult to compute 
in most cases and the weighting factor c is also difficult to 
Calculate. Moreover most of the work in univariate random 
number generation has been aimed at avoiding having to calcu- 
late inverse cumulative distribution functions such as 

Be (-) and eat) These are the reasons why many proposed 


methods are specific to a specified marginal distribution. 


Again special properties of certain random variables such 

as infinite divisibility have been exploited to give easily 
generated bivariate random variables, often though with 
limited ranges of dependency. One very clever scheme by 

Gaver (1972) to generate bivariate exponential random varia- 
bles uses the fact that the sum of a geometrically distributed 
number of exponential random variables (Y) is exponentially 
distributed and that the minimum of this geometrically dis- 
tributed number of independent logistic random variables (Z) 
is exponential. Clearly when Y is large, Z is small. This 
scheme is of course very specific to exponential marginal 
distributions and, via an exponential transformation, to 
uniform random variables. To avoid these kinds of limitations 
and to make the generation of bivariate random variables 
Simpler and more automatic in simulations we develop here a 
scheme presented in Jacobs and Lewis (1977). 

This scheme, the mixture-truncation method, is a very 
general tool which requires only that a method be available 
for generating random variables with the desired marginal 
distribution. The mixture-truncation method scheme for 
generating bivariate random variables is as follows. 

mete (x) be the common marginal distribution, of the 
bivariate random variable (Y,Z). Let o be the desired 
correlation between Y and Z (which may or may not be attain- 
able generally or in particular with the mixture-truncation 


scheme). 





Define the transition matrix P as 


Oy L-a, 
4 = 
l-a, O 6 
Wieiestansonary vector 
oe = 1l-a, l-a, 7 
= = L-a,t+1l-a, l-a,+1l-a, 


and let the range of the random variable X with distribution 
function F(x) be x. Then generate (Y,Z) as follows. 
ieee (initiallization) 


1) Choose an "allowable" Xo £rom (x,,X,) (the 


allowable range of Xo will usually be smaller 
than x and depends on o and F(x)). 


Ha) Set ™) = F(x), 1) = l- 1, and compute Oy and 


& 


3° 
iii) Denote by xX) the random variable X truncated to the 


é 


left of x, and by X, truncated to the Pag hiieod 


Zz 


Xe 
O 


2. (Generate Y¥). Choose Y from X, with probability Ty) 


i 
or choose Y from Xo Wath oropabl lity Ty 

3. (Generate Z) 

1’ choose Z@ from x) with 

probability a, Or choose Z from Xo with probability 


je Lf Y is chosen from X 


1-4. 
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MP heer ats Chosen From X51 choose Z from xX, with 


PBobabi lity l-a, or choose Z from Xo with 


probability Oo 


In Section III, we will show that the correlation » between 


Meand 2, if X, is properly chosen, is 


» = ee] a ag 
= cee, 
where 
B = a, - (l-a,), 


(u, - 45)? 
eee 1 2 
\4=- p-OQOOs>"  _,, 
2 
a 


= E{X,}, 2D = E{X,} t 


Ox = VAR[X] , 


Moreover, the generated bivariate random vector (Y, Z) 


wale 


have marginal distributions F(x) and correlation po and the 


feentee distribution function of (Y,Z) will be 


Pe 





et 
— F(z) F(y) peace, Y < xX 
ae 


(l-a,) 
fo, + 7 [F(2) -F(x,) J} F Cy) Te caer YX 
EB Yiiz) Meo} 
tay + SE CY) -F (0) 1 (2) LE weeks Y > X 


T 
Hl ; 
a m, + (1-a,) [F(2)-F(x,)] 7 pa Z<Xo, yx 


r 


i 


c ((1eag) + Z2OF (2) -F (xg) 1} OF (Y) -F (xg) 
These relationships will be developed in Section III. Note 
Eiga t. (Y,2Z) may be continuous or discrete random variables 

or mixtures of both, though in this thesis we concentrate on 
continuous cases. 

The key problem in this very simple algorithm comes at 
the initialization steps i) and ii). There are two degrees 
@eeeeedom in the selection matrix P but setting ™ = F(x) 
constraints reduce to one degree of freedom. Specifying 
a desired correlation further constrains the degree of free- 
dom, though not completely. Subject to the constraint that 
ay and a5, are probabilities there may be 

a. no values Se which will give (Y,Z) with correlation 

p 


ey one value Xo 


Cc. a range Xo" 


me 





The main numerical problem of initialization step (i) then 
is to compute the range of allowable x for a given p and 
B(x). 

The main statistical problem is then to choose which of 


the bivariate random vectors (Y,Z) indexed by x_ « ee use. 


O 
An alternate solution to the statistical problem, giving 
another algorithm is then to let XS have some distribution in 
the range x, possibly uniform or triangular, this not only 
alleviates the problem of picking a particular x, « xe but 
it also smoothes out the distribution of (Y,2) and possibly 
makes it continuous. The computation of x, for given p is 
illustrated in subsequent sections for uniform, exponential 
and gamma marginals for (Y,Z). 

Before doing this and developing, in Section III, the 
results already given here, we review a few existing 
methods for generating bivariate random variables in Section 
II. These are methods which seem fairly tractable, for use 


Pama COMpPuter. 
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fie REVIEW OF EXISTING METHODS 


There are several existing methods for generating 
bivariate random vectors, but most of them are specific to 
Sarmiwemlar Marginal distribution and use inverse transfor- 
mation method. The inverse transformation method is a very 
useful univariate procedure which, unfortunately, is not 
possible to use with many distributions because it is diffi- 
cult and/or uneconomic to compute the inverse functions. We 
survey here some of the methods which are germane to this 
Baesis, in particular concentrating on bivariate random 
variables with uniform, exponential and gamma marginals. First 
we review the problem of determining the range of correlation 
coefficient p which can be obtained for bivariate distribution 
with given marginal distributions, again considering only 


the continuous case. 


A. RANGE OF CORRELATION COEFFICIENT op 

suppose that Y, Z are random variables with an arbitrary 
jJOint distribution F(y,z) with finite second moments. Then 
in general the correlation coefficient p can take all values 
in the closed interval [-1,1]. But with specific marginal 
aor bution F, Cy) and Fo(z), the class of all F(y,z) need 
not attain the values of -l, 1, of p. The necessary and 
sufficient conditions that there exist determination of F(y,z) 


with p equal to 1 and -1l are given by Moran (1967) as follows. 


14 





mmetchnere exXiStmconseants @ and 6 such that aY+8 has 
the same distribution as Z, 
ieee the distribution of Z is symmetrical about its 
mean. 
To see this, rescale Y to have the same mean and variance 
as Z. If there exists an F(y,z) such that p = 1 we have 
E(z-y]° = 0, so that Y = Z with probability one. Then Y 
must have the same distribution as Z. On the other hand 
if there exists an F(y,z) such that p = -1l we shall have 
B[Z+Y¥]° = 0, and Y = -Z with probability one. Thus if both 
bounds are attainable Z must have a symmetric distribution. 
Given general marginal distributions Fy Cy) and F(z), 


Mardia (1970) showed Frechet bounds as 


max[0,F, (y)+F,(z)-1] <1 (GAS minlF) (y),F, (2) ] (II-A-1) 


From this we can find the range of possible values of op. 

For simplicity we now confine our consideration to distribu- 
tions F(y,z) of positive random variables whose deriva- 

tives Fs Cy) and Fs (z) nomcinianGtiy = DOsdtaive for y > 0, z > 0, 
respectively. Suppose also that the variates are scaled to 
have unit variances. Let G, (u), G, (v) be the inverse func- 


tions of Fi Cy), Fo(z), i.e. 


F,(G, (w) ] = WwW 


where 0 ieee 1, 2. Then the correlation coefficient 


in 





between Y and Z is given by the equation 


oO oO 


fy z aF(y,z) - E[Y] E[Z] (II-A-2a) 


OD 
lI 
we, 


0 0 


Tee 
J ,! G, (u) Gj(v) aK(u,v) - EY] E[Z] (II-A-2b) 


0 


where K(u,v) is the joint distribution of the quantities 

U = Fi ty), V= F.(z). Then U,V are jointly distributed on 
meemogquace 0 < U< l, 0 < V < l, in such a way that the 
marginal distributions are uniform on the unit intervals. 
From expression (II-A-1) the minimum correlation is attained 
when the probability is concentrated uniformly on the line 


U+Ve= 1. The minimum value of po then is 


i 
eel = e G, {u) G,(1-u) du - Efz] Ely] (II-A-3) 


The corresponding F(y,z) will be a singular distribution with 
all the probability concentrated on the line Fi(y) + FP. (z) = l. 
In fact Y and Z are an antithetic pair. The maximum value 

Of 9 is attained when the probability is concentrated uni- 


formly on the line U = V. Then the maximum value of p is 


1 
A = ») Secs edie e 2) IY), (II-A-4) 


Et 





with the corresponding propability concentrated on the line 
F, (y) = F(z). 

By using this result we will get lower and upper bounds 
of the correlation coefficient for uniform, exponential and 
gamma marginal cases. In the uniform marginal distribution 
case, we can see o0_. = -l and o = 1 from Moran's condition 

min maxX 

ii), i.e., uniform distribution is symmetric about its mean. 
For the exponential marginal distribution case suppose Y and 


Z have exponential distributions with unit mean and variance. 


Then 


= Beers ay, 2 ~-” 
Fy Cy) = | e-, F(z) al! e 
and the inverse functions are 
G, (u) ernie — il), 
G, (v) = - Iin(l-v). 
From equations (II-A-3) and (II-A-4) 
a 
oie ) G, (x) G,(1-x) dx - E(¥] E(Z] 
1 2 
= fin x 1n (1l-x) dx-1 = 1-5 
0 6 
> -0.64493 


av) 





1 
J G(x) G(x) dx - E[Y] E[Z] 
0 


O 
ll 


Max 


i 
feos in x dx —- 1 = 1 
0 


The oe can be attained when all Y and Z are concentrated 


on the line e” +e% =1. Also the Onax Can be attained 


= 2 
’ 


when all Y and Z are concentrated on the line e” =e 
1.e., Y= 2. 
For the gamma marginal distribution case, suppose that 


Y and Z have gamma type distributions with probability density 








functions 
f, fy) = = e Y iia CeOG my 0 
poe ay* 2 %@>0 ,220 
Then 
BY) = a, E{Z] = 8, VAR[Y] = a VAR[(Z)] = 8 


The oe . Will be attained when the probability density is 


concentrated on the line 


l Y -x a-l l -x 3-1 





HS 





This defines y uniquely as a function of z which can be 
written y = A(z). The ona is then, on rescaling the 


covariance 


eo 


onin = u A(u) £,(u) du - o8}/(a8) 


nee 


When a,8 become large, op... tends to -l. And the o 
min max 


will be attained when the probability density 1s concentrated 


on the line 


1 Y -x a-l a -x a-l 
Tey = x ax = J e x ax 





This defines y as a function of z which can be written 


y = B(z). The oer ee then is, on rescaling the covariance 
Q = { fu Be) f(y) see — a8 }/ (ap) 2’? 
max 0 1 

when a = 8, Oat tends to 1. Schmeiser and Ram Lal (1979) 


Showed the obtainable correlations between random variables 


Y and Z having gamma marginal distribution with density function 


a,7l 
f(x) = L(x/8 5) exp(-x/8;)/(8;T(a;)] 


mex > QO, a; > 0; Bs POs. OL = Ly 2% 


JUS) 





The Figures (II-a) and (II-b) show the obtainable 


Perrcelation as a function of Oo given a, = ances eho. 





OS) 


Figure II-a. Obtainable correlation as a function of 


o5 Lor a, = ae By = 1, by Schmeiser 
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Beh oe a 0 Ve ha ee eee rp 
. a ‘7 eae Pk ta eh « : 
Te n ike tank, 


wa te’ 
wes 
er 4 


* 2 


p ~ 


ee ee ee on ae Ge ee on ae on on an on on an om on Ste 





Figure II-b. Obtainable correlation as a function of a 


e@ te a, = 5: 84 = 1,by Schmeiser 


Z 


Ze 





B. GENERAL METHOD REVIEW 
1. Johnson and Tenenbein's General Method 
A general method of constructing a bivariate dis- 
tribution, whose marginal distribution functions are B(x) 
and Poly), 1s proposed by Nataf (1962), can be represented 
as follows. 
General Method 
1) Consider any two continuous random variables U 
anewy Wicehepmebability density functin h(u,v). 
ii) Let X' = H, (u) ana Y* = H,(v), where H, (u) and 
H. (v) are the cumulative distribution functions 
of U and V, respectively. 


1ii) Define 


E: wil ' = ek ae 
ee eo ee ela) | (II-B-1) 
and 
Y = Fit (y') ee) (II-B-2a) 
2 DD 
Or 
ee =e ie yy) (z= B= 2b) 
2 2 2 


2 





Since X', Y' and 1-Y' are uniformly distributed over the 
interval [0,1], X defined by expression (II-B-1l) and Y 
defined either by (II-B-2a) or (II-B-2b) will have a joint 
distribution whose marginal distribution functions are Fy (x) 
and F,(y)- The method is probably what one would think of 
first but again involves inverses. Based on this general 
method, Johnson and Tenenbein (1979) developed two proce- 
dures for generating (and also constructing) bivariate dis- 
tributions whose marginal distributions and measures of 
dependence, as given by Kendall's F or the grade correlation 
coefficient Pg, are specified. Note that these Kendall's fF 
and grade correlation coefficient po, are not the same as the 
Ordinary product moment correlation coefficient p we have 
been considering; these measures are discussed by Kendall 
(1962) and Kruskal (1958) and can be defined as follows. Let 
X and Y be continuous random variables having some joint 
Bropability density function. Let (X),¥,), (X,Y) and 
(X5,Y3) be three independent pairs of observations having 


the same joint density function, then 
aed 4 —X.) GY =Yo) > 0) = 1 
p =) 6 P[(X,-X5) (¥ aaa) OS 


S ics 


The first procedure for generating bivariate pairs, called 


the WLC (weighted linear combination), defines 


BE 





ce Oe Cli-B- 2) 


and 


Vee eeu. + (=e) y- Ce 5-4) 


Seeemor- c < 1, U' and V' are independent and identically 
distributed random variables with probability density func- 
tion g(-). Then X and Y are obtained from the general 
method using equation (II-B-l) and (II-B-2a) if the depen- 
dence measure is positive or equations (II-B-1l) and (II-B-2b) 
if the dependence measure 1S negative. 

In order to apply the WLC procedure, we must obtain 
expressions for H,(u), H(v), 0 (u,v) ancduu cue) pe ineaterms 
of c and g(-). The values of H, (u) and H. (v) allow us to 
apply the general method for a given choice of c and g(t). 
The expressions for T(u,v) and 0. (u,v) allow us to specify 
c for a given choice of g(-) in terms of the required value 


of either T or O° From equations (II-B-3) and (II-B-4) 


u 
H, (a) = - giceece (II-B-5) 
apes i) GCA) gigal) dusdy" ~~ (I1-B=-6) 
R 
1 


24 


where 


Ry CC )me el. + (l—c)y" < y} 


amcaethe jOint density function of u and v is 








h(u,v) = peo glu) g (PS) (II-B-7) 
mev) = 12 f __f Hy (ai, (w)h (u,v) duay - 3 (II-B-8) 
fu) «= «| eee ec eee (II-B-9) 
where 
g(t) = _f st + t) g(w) dw 
eGo _f 9g dx 


Using equations (II-B-5), (II-B-6), (II-B-7), (II-B-8) and 
(II-B-9) , we have to evaluate H, (u), Ha (v), Pinel i tela) a, 
and p, (u-v) for all values of c and for specified g(°:). 

As Johnson and Tenenbein mentioned in their report, these 
integrals are quite tedious to perform, so they gave some 


computation results in their report. 
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The second procedure, called TVR (trivariate reduction) 
is discussed by Mardia (1970). In this case U and V are 


defined as 
ee — ee ees Ao (II-B-10) 
Vo set | (II-B-11) 


memo < ©, and U', V', and Z' are independent and iden- 
tically distributed random variables with probability density 
Mounction g(°). 

Like the WLC procedure, one needs to define H,(u), 
H(v), Reta, V) ; p (u,v) and 7T(u,v) as a function of 8g and g(°:). 


From equations (II-B-10) and (II-B-11) it follows that 


H,(u) = H(v) = J Jaa") gz) du! dz (II-B-12) 
3 
where 
R, = i(u',z'): ul + Bz' < ul, 
Mu, Vv) = i g(u- 8z)g(v-86z) g(z) dz (II-B-13) 
and 


Z6 





co 


mee sie (ee) I g(t) dt — 1 for Bl) 


where 


co 


{ g(w+t) g(w) dw 


-_ © 


gq (t) 


oO 


fF G5(x) dx 


G. (t) 


using equations (II-B-12), (II-B-13), (II-B-14) and (II-B-8) 
the same computations are needed. 

In this general WLC and TVR methods, there are two 
problems. One is the need for inverse transformations. The 
other is the need to generate coupled random variables to 
create dependence after the inverse transformation is applied. 
Doing this by linear combinations is not necessarily the 
Simplest and most felicitous with respect to calculation of 
correlations and range of correlations. 

A significant simplification can be achieved by 
obtaining a "smooth" bivariate uniform pair (U,V) whose 
correlation can be varied between the limit of correlation 
which can be obtained for uniform marginals. These limits 
are (-1,1) since the antithetic pair (U,1-Y) has correlation 


=a 
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Note that the problem is essentially one of obtaining 
: ; ome : 
' a bivariate uniform random variable with positive correlation 


because 


= -V 
Py 1-V COrr (h,1—-V} 


corr(U,V) = ig ae 
a 


Seu, v = corr(1-U,V) 
Two fairly simple schemes for obtaining bivariate uniform 


mrstcribucions are discussed in the next section. 


Zeer evilEW OF SOME SPECIFIC METHODS 

The mixture-truncation method is a general algorithm 
for generating bivariate pairs (Y,Z) with given marginal 
distributions and requires only the availability of a generator 
for the marginal distribution and calculation of constants. 
Unless randomization is used, however, it does give a bi- 
variate distribution with a discontinuous joint distribution, 
as in Section (III-D). This could be a disadvantage. It 
1s interesting therefore to compare it to other generations 
which are specific to the given marginal distribution and 
we do this here for the cases for which the mixture-truncation 
method is illustrated in this thesis, namely uniform, exponen- 
tial and gamma. 

There are many schemes available for uniform and exponen- 
tial cases though it should be noted that smooth, simply 
generated pairs with easily calculated correlations have 


only recently been available. In particular we give here two 


Ze 





very recent schemes for smooth bivariate uniforms which are 
easily generated and whose correlations can be calculated; 
similarly two schemes for exponential which are essentially 
the uniform schemes after transformation. When one comes 
to gamma marginals one is again in difficulty and we describe 
a recent scheme by Schmeiser and Ram Lal which involves 
however, computation of the inverse gamma distribution 
function and very difficult initialization computations. 
The mixture-truncation method is quite competitive here. [In 
fact if one can compute the inverse gamma distribution then 
bivariate gamma's can be computed as [Y = oe Z = ae 
where U and V are a bivariate uniform generated by the schemes 
mentioned above; this would be simpler for any marginal dis- 
tributions than some of the general schemes mentioned in the 
previous sections. Comparisons of the generating schemes given 
in this section with the mixture-truncation method are given 
in Section IV-D and Section V-D. 
1. Lawrance and Lewis's Bivariate Uniform 

Lawrance and Lewis (1979) show that a simple trans- 
formation of the NEAR(1) (New Exponential Autoregressive) 
process gives a two-parameter family of Markovian random 
variables with uniform marginal distributions. This method 
generates acorrelated uniform pair as a multiplicative mix- 
ture of uniform variables, where the parameters a and g take 
on values between zero and one, with the condition that they 
not both be one. Let Y be a uniform [0,1] random variable, 


and define 


2g 





ee WD a 
vf = 

E Se 

where 
U Wp ae 
1 - (l-a) 8 

Ee — 

aloe aB 

S ee acne 


where U is an independent identically distributed uniform 
(0,1) random variable. The correlation structure is defined 


as follows as a function of a and 8; 


3 a B 


0 Sora Ts ea 
Y,2Z eee oe bee Coc 


The model can be reduced to a one-parameter model by suitably 
relating a and 8, e.g., a = 8 and all positive values of 0 
can be obtained. Consequently Y and 1-2 will have a full 
range of negative correlations. A generating procedure for 


this bivariate pair of uniform random variables is as follows. 


Generating Procedure 
le Mlnitialization). For given correlation, find suitable 


parameter values a and 8. 


as 


set Y < U eran 


2. Generate Ul, li r= l-a, P < 


3. Generate U.- 
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U 
Peete: U, < ©, set U + = and go to 8. 





2 
U,-2r 
5. Otherwise U < eae: 
B U 
6. If U< P, set 42+ Y° Jap and go to 10. 
7. Otherwise, set Z < ye (Fae) FF and go-to 10. 


Pa: UU < P, set Z + = elavel (ee) ces) I) 


eee 
1-P . 


9. Otherwise, set Z < ( 
10. Go to 2 and repeat until the desired number of bi- 
variate uniform variables (Y,Z) are obtained. 

The program is listed in the appendix and the scatter plot 
and bivariate histogram of resulting output are in Section 


mV-D. 


2. Bivariate Uniform By Transformation of Gaver's 


Bivariate Exponential 


In Gaver's bivariate exponential to be discussed in 


Section II-c-4, we can define a bivariate exponential random 


vector, (Y,Z) with correlation p = - P/2 as follows: 
i Pp 
y= 1in(——— - -—] (II-C-1) 
yl/N (y_p) 1=-P 


and Z is, conditional upon N = n, a gamma random variable 

with shape parameter n and mean n(l-p), where N is a geometric 
random variables with probability density function 

f(x) = mil-p)*-, S-ele 2, es. and U 1S a uniform random 


Variable over [0,1]. From this we can generate a pair of 


Sol 





uniform random variables (V,W) by transformation. Since we 
know that to generate an exponential random variable X from 
a uniform random variable U[0,1), we use the following 


transformation: 


ou] s= Int le— U) (II-C-2) 


The resulting X has the exponential distribution with unit 


mean. From equations (II-C-1) and (II-~C-2), we can generate 


V, a uniform random variable over [0,1]: 


pas i) yl /N 
1- pvui/N 


Further also we can generate W, a uniform random variable 


over [0,1]: 


This follows since we know that Z has, conditionally upon 
N=n, the gamma distribution with shape parameter n and 
mean is n(l-p). 

In general the resulting V and W will be negatively 
correlated and (V,1-W) will be a positively correlated pair. 
Because the correlation structure will be changed by trans- 


formation, the resulting (V,W) need not have the same 
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correlation with (Y,Z). The correlation coefficient of 


(V,W) will be computed as a function of P as follows. 


o0 = 12 E[vw] - 3 


where 


E [Vw] 


i 


E, (E{VWIN=n]] 


iE - 17a aie 
= EI j {tS ee tal ae du}"] 
0 


ee 0 
Unfortunately this computation of o as a function of P is 
difficult. As an example we used here the same function for 


p as holds in the exponential case. 


Generating Procedure 
meee Lnitialization) 
a: ) Compute P for given correlation o 
ii) Choose N from a geometric distribution with 


parameter 1-P 


2. Generate a uniform [0,1] random variable Uy and define 


ani 
1 


_ IL ie 
1 |S UL 


(=P). U 


3. Generate N = n uniform [0,1] random variables U5, 


jme- 1, ..., nN, and define 


a 





fmmebeliver (V,W) and go to 2 until a sufficient number 
of bivariate pairs have been generated. 
The program is listed in the appendix and the scatter plot 
and bivariate histogram of resulting random vectors are in 


Section IV-D. 


feeevarshall and Olkin'’s Positively Correlated 


Bivariate Exponential 
Suppose X is the age, or length of service of the 
first device at the time of death, governed by two indepen- 


(t) with parameters \ 


dent Poisson processes Z,(t), Ll’ 


ane 


A419 respectively and Y the same of the second device, 


governed by two independent Poisson processes Z(t), 219 (t) 


with parameters \ respectively. Further suppose the 


ee 
second device is placed in operating after a time 6 later 
mien the first device. Then the joint distribution of (X,Y) 


1s defined as 


Pane > X,Y > y) SEG) 


exp{-\)X-A oy — A, 5max[x,ytmin(x,6) ]} Li oe 


exp{-\,x~) max{y,x+tmin(y,6)]} Deen VG 


J A 
Cc) 


ee > 


Cini C=.3) 


If 6 = 0, then equation (Ii-C-3) reduces to 
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mie x, oy) = F(x,y) 


= exp{-A)x - oY - AWo Mmeaccts< 7 )) 5 


Note that since the marginal distributions do not depend on 


Oy 


i 
E[X] = ’ 
eet A312 
(A, Ayo) 
l 
Ea 
No FAd9 
VAR[Y] = : 5 
(A5 + Aqo) 
We also have, for 6 > 0 
. ee 
l 12 


B[XY] cn a Pn 
Sem oe yA FA) 45 4D) 


and the correlation 
r ~-(A,+A..)6 


eomr(X,Y) = o = —e t 12 (II-C-4a) 


where 


chs 





= 

rh 
Or 
it 


O, then the equation (II-C-4a) reduces to 


0 = —— (II-C-4b) 


This characterization can be represented in terms of 


mi@eeendent random variables. For 5 > 0, if there exist 


independent exponential random variables Uj, Uo, U, and U, 
with respective parameters hye hos Ayo! Ayo: then 
xX = min(U,,U,) 
min(U.,U,-5) af U, Zt 
Y = 
min(U,,U,) at U, ane 


It can be verified formally from the relation 


P([X>x,Y>y] = P[U, >x,U,>yI{P[U,>x,U > yt6|U, > 6] 


3 


eP[U,>6]} + P[u, >y]P[U,>x|U, oe la 


In the 6 = 0 case, the representation yields 
xX = min (U,,U,) 
Y = min(U,,U,) 
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This bivariate distribution has a line-discontinuity, if 
E[X] = E[Y], along the line x = y, 1.e., in the case where 

x =v U,- Thus this bivariate exponential is not as smooth 
as others. In particular the NEAR(1) process of Lawrance 
and Lewis generates bivariate exponentials with a continuous 
mememty cunction. But if we consider the case 6 > Q, the 
line-discontinuity will be removed. In either case, the 
resulting (X,Y) pairs always have positive correlation as 
shown in equations (II-C-4a) and (II-C-4b). If we use these 
methods to generate positively correlated random vectors 
(X,Y) for exponential marginal distribution with unit mean 

A, and A LOnsgiven 


ie <2 a2 


correlations to generate independent exponential random 


and given op, we have to compute A 


variables. For simplicity we will consider the method with 


§ = 0. From the relationship A = Misen Noe) 


2 ee 


1 
E{X] = = J] 
Ay A192 
a 
E{Y] = = |] 
ho a 12 
eee 
O X ’ 
we can compute hye ho and di0 as functions of o. 
Mae = Zay 1 aa) 
Bape oO AT 
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Generating Procedure 
fee (initialization). Compute Aye r, and 12 for given 
Oop < tL. 
2. Generate three independent exponential random variables 


Eeawanad E. with unit mean. 


ae > 2 


3. Rescale the independent random variables and set as 


U, + ES, 


U, <a Eo/\5 


U3; + E3/h1>5 


4. Define X and Y as 


a min[U, ,U] 


K 
ll 


min[U,,U,] 


59. Deliver (X,Y) and go to 2 until a sufficient number 
of bivariate pairs have been generated. 
The program is listed in the appendix and the scatter plot 
and bivariate histogram of resulting random vectors are 


Shown in Section V-D. 
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eecaver os Negatively Correlated Bivariate Exponential 


This negatively correlated bivariate exponential 
is generated from the following considerations. Suppose a 
particular system has N defective elements; here we assume 
N has geometric distribution with probability density 


munction 
x-] 
foc me = 4) le, 


The time to failure of each defective element is T. with 
density function F(t). Further suppose repair is carried 
out after all defective elements are discovered. If we 
Meetine R as repairing time which is distributed as nce 
eoovided N =n, and define T as minimum of T5 then the 


Beeeiesdistribution function of T and R is obtainable from 


co 


Pie t,Rsx] = ) (1-F(t)) PP. G (x) 
i ne) 
where 
n-1l 
om = (l-p) p pe ae. <a elie 2 | 
Then 
E(R|T>t) = >is 


i= Pipe FC) | 


3h) 





where E[{R] is the expected value of an element repair time. 


Moreover, 
B(R/T=t] = E(R] jee 
= E(R] (1 + p5(1- 9(t))1 
where 
ae) = 6 PIT<tj 6= F(t) 


lee be 


Brom this regression, we know that if T is long, R is short 
and vice versa; the negative relationship is clearly present. 


Purther we can calculate Cov[R,T) as 


Eon eae 


deve 2 ENR) EIT) = 


nO, 
where T(1) has the distribution of the smallest of a sample 
meetwo from ¢(t). 

It will not be shown that one can select F(t) = 1- F(t) 
in such a way as to induce exponentially distributed time 
Memeystem failure, T. Since we know that exponential G, 
geometrically compounded, yields exponential R, the outcome 
will be a (T,R) pair with exponential marginals and negative 


correlation. 
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Coe! —- Fe) 
satisfies 
. = n on-l -At 
ence Means (l-p) = e 
n=l 
mien £f is exponential with mean 1/dA. Solution for F 
yields, 
= _ 1 
F(t) = —— ed 


p+ (l-p)e 


which is a logistic distribution, left truncated at t = 0. 


meee — 1, then T is exponential with unit mean, while 
E(T ii)! = 1/2, so 
corr(R,T) = - & E[R] EIT] 


If G is chosen so as to make R exponential then it may be 


Shown that 


IL 


PS osc a] 
2 E[(N] 


COEE(R;T) = - xll ~ 
Consequently, a greatest lower bound for the correlation in 


this model is - 5. 
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Generating Procedure 
ee ieee 0 < S055, Set P <« -20 


2. Generate a geometric random variable N with parameter 1-P 


3. Generate a uniform (0,1) random variable U and define 


Y as 


4. Generate a gamma random variable, Z, with shape 


parameter n, and mean is n(1l-p) 


Seeeebeliver (Y,Z) and go to 2 until a sufficient number 
of pairs are obtained. 
For obtaining negatively correlated bivariate exponential 
random vectors, this algorithm is very simple and one of the 
Few available. Moreover the correlation is known explicitly. 
The program is listed in the appendix and scatter plot and 
bivariate histogram of resulting random vectors are shown 
in Section V-D. 
5. Arnold's Vibariate Gamma Generator 
Arnold (1967) developed a Trivariate reduction method 

to generate bivariate gamma random vectors having positive 
correlation with 


COLE >. < min(a,,05)/(a,a5) 7/7 ; 
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where O45 and a, are the given shape parameters for the marginal 


2 
distributions. Letting "gamma (a,8)" denote the gamma dis- 
tribution with shape parameter a and scale parameter 8, 

so that 


_ } poecus * 
f(x) = ar tay (3) exp (-x/8), 


where [(a) is the gamma function and 
E(x] = a8 j VAR[x] = ag , 
the trivariate reduction algorithm proceeds as follows for 


given shape parameter values A,r A and 9 to generate unit 


scale gamma variate i.e., 8 = l 


Generating Procedure 


2 
1. Generate Gy ~ gamma (a, - 0 (a,a5)*/ aa) 
72 
2. Generate G, ~ gamma (a, - 9 (a, a>) 4 vals) 
3. Generate G, ~ gamma (0 (a,a)°/*,2) 


4. Define Y and Z 
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These Y and Z are unit gamma variates with 8 = 1, but Y and 
2 can be multiplied by By and Bo, respectively, to obtain 
any desired scaling. In the algorithm, G4 is a common com- 
ponent of both Y and Z, thus inducing positive correlation. 

When a positive, but not extreme, correlation is 
needed, trivariate reduction is an excellent algorithm, 
using univariate gamma generators. If the marginal distri- 
butions have the same shape parameter value a, then the 
SeeretatiOn limitation is 0 < p < 1. But in the other case 
the upper limit is smaller than l. 

Note that this algorithm exploits the infinite 
divisibility of the marginal gamma variables and is applicable 
to any infinitely divisible marginal. It is commonly used, 
for instance, to get bivariate Poisson random variables. 

6. Schmeiser's Bivariate Gamma Generator 

Schmeiser (1979) developed a family of algorithms, 
any member of which can generate bivariate gamma random 
1" %» and allowable 


correlation p. In this section we will discuss his general 


vectors having any shape parameters a 


procedure. 


Be Z, Wy and Wo are independent gamma random varia- 


bles with shape parameters r, and 6 respectively, and 


S47 De 
U is an independent uniform [0,1] random variable, and either 


V= U or V = 1-U; then 


We F-*(U) +2¢W eae) 


1 I 
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1 and FY are the inverse distribution functions 


where F 
Se gamma (n,l) and gamma (m,l), respectively, distribution 
function. The resulting Aas X. are gamma random variables, 
with shape parameters 


a ee Obie 
(II-C-6) 

O = mtr et S51 

respectively. This result follows immediately from the 

reproducibility property of the gamma distribution and from 

Moting that F ~(u) and Peo (leu) are each random variables 

having CDF F. This scheme generalizes the trivariate 

scheme but brings in the inverse CDF. 

The correlation coefficient is defined as 
B(F + (u) Pol (y)] -nmte 
= n m 

Py x i \75° °° iii se="7 ) 

ee? (a, a5) 

The remaining problem is to select values of the five 

Parameters n, m, r, 6, and 6 


JL 


distributions and correlation. The following conditions 


5 to obtain the desired marginal 


must be satisfied. 


45 





me E+ 6 


i 
Q 


2 2 
E(F + (u) Ft (y)] - (nm) + rx 
n m fy 
aa Bag = =k (II-C-8) 
(aa 
ee 2 
His Ho, Te Say 65 oan 


Since we are uSing five variables to satisfy three equality 
conditions, finding a set of parameter values corresponds 

to finding a feasible solution, rather than an optimal solu- 
tion, to a nonlinear programming problem. 

An efficient solution procedure for determining 
parameter values is important, since substantial computation 
is required to determine whether or not conditions (II-C-8) 
are satisfied for given parameter values. Most of the 


computation is involved in calculating 


a E[F>*(u) Fo*(v)] -nm-r 
(a,a,) 74 q 


Since the expected value must be calculated numerically 


using any one of the following three integrals: 


aaa -1 
eee) Eos (a) du (II-C-9a) 


ce 


f FOO (BP (x)) x” exp(-x) dx/r(n) (II-C-9b) 
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co 


f Fo*(F_(-1n y)) (-1n y)® ay/r(n) Gai) 
0 


1 ew) en 


jeep > O. If p < 0, then replace Fo*(u) with Be 
equation (II-C-9a), replace F(x) with 1- F(x) in equation 
(II-C-9b), and replace Fo (-in y) with 1-F,(-In i) baeeee alt lon 
mee-C-9Cc). 

Schmeiser seemed to have best results in terms of 
a subjective trade off between speed and accuracy, using a 
24 point Gaussian method for the integration. He selected 
the parameter values from the feasible values satisfying 


conditions (II-C-8) by making the curves of regression 


E[X,|X,] and E[X,|X,] behave as desired. 


Generating Procedure 
1. Generate x) ~ gamma (n,1) 
2. U~<-+ PS fXy) 
Miro < 0, U<+ l-u 


3. Generate Z ~ gamma (r,1l) 


4. Generate Wy ~ gamma (6,1) 


5. Generate Wo ~ gamma (65,1) 


6. XX, + ay + Z + Wy 
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Based on these procedures he developed a family of algorithms 
which can provide variates having any theoretically possible 
correlation op. Anyway, for gamma marginal distributions, 
not all correlations are consistent with particular shape 


parameter values. Schmeiser shows the obtainable correlations 


mea function of Os given a, = 1 and 5 as shown in Section 
Tf-A. Note that only when 1 = o5 1s it possible to obtain 
0 = l- Likewise p = -l is not possible except in the limit 


as a, and a. tend to infinity. The maximum and minimum 


Z 


possible correlations, given in Moran (1969), occur when 


-1 
X =e et CX- )) 
2 a5 Oy a 
and 
ee (LF (X,)), 
Ae a 


respectively, where PX) and P+ (1) are the cumulative 
distribution function and inverse CDF, respectively, of the 


gamma distribution with shape parameter a. 
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III. GENERAL MIXTURE-TRUNCATION METHOD 


Denote by (Y,Z) the bivariate random pair, where each 
has identical marginal continuous distribution F(x), and 
denote a general random variable from this distribution by 
X. The argument is not specific to continuous random varia- 
bles; this aspect comes in only in the computation of the 
correlations and can be developed in a parallel fashion for 


discrete marginal distributions. Let 


(T1775) ’ 


and let be an X truncated to the left of a fixed point Xos 


x i 
Xo be an X truncated to the right of X51 SO that 
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- = PIX, < x] = 
F(x) - F(x) | 
Ge): Te aoe aS x 
Oo 
In addition we set nr, = F(x_), t~ = L-7wr, and choose 
it O 2 JL 


* and Z as follows. 
1) Choose y from Xy Wathen pie @bals 1 lartey. Ty and then choose 


or from x%. with 


7. Lrom Xy Wolick?, jsve@leysi aii Ih sie Oss D 


Prepability 1l-a,- 
11) Choose Y from X., Wate or Obab il Tey Tos 


and then choose Z from xy with probability “l= «a 


where +17, = l, 


9! or 


from x, with probability a 


5° 
If we choose (Y,Z) as in the above procedure, then we can 


make the following two theorems. 


A. MARGINAL DISTRIBUTIONS 
Theorem l. 
The marginal distribution of (Y,Z) becomes F(x) for both 


mand Z. 
Proof 


Peeemarginal distribution of Y 


By definition Y is the mixture of X, and X, with proba- 


Ih 2 
be lity Ti To respectively. That is 


Ty) F.. (x) + T 5 FY (x) 


ee (3c) 
: 1 2 
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ae F(X) 2 Oo 
EF — 
y **) 
F(x) - F(x) 
TT - l + T I+ ¥FTx,)_ Peo xX. > x, 
But since we define a F(X,)3 tT = 1-1, = 1- F(x), we 
have 
F(Xx)) 
TT = ob) io << k 
le ig — “Oo 
1 
Py, (x) = 
F(x) - 7) 
TT) + T, = F(x) i xen OX 
2 
= ec) in all cases. 
So Y has the marginal distribution F(x). 
feeeMarginal distribution of Z 
imtiey is from Xj then 
Seo a. FY (x) + (l-—ayz) FL (x) 
aT i Xy x5 
ees ees) ° 6 if x< x 
i F(x.) 1 ae) 
F(x) - F(x) | 
ol cian E(x.) if x > xX, 


Sl 





mM Y is from X then 








De 
mtx) 6COU=)3=S s (1 —~ a) FO. (x) + aw F. (x) 
Z. Z x) Z Xo 
cS 
Bex : 
(l- a.) F(x.) - 5 0 ve x 
F(x) - F(x_) | 
(l- a.) 2 Lose a  L- F(x.) ae <a ao 
50, 
eC = CUCU, OF OC XY) hUv+thCUt. 6FQ OCC x) 
Z ik 2 Z Z. 
TO) a + tm, (1- a.) ao Eee XS 
O O 
Pei eee | 
T, (a, + (l-a,) eee ae) Le 3 co 


F(x) - F(x)) 


oo =a 
O 


+ tma((1- a.) 


and we defined ™) and T > as follows. 


il l-a,+1l-a, 
T = oe 
2 l-a,+tl-a, 


a2 





From this, we know 


eo th e)) 
2 
If we use this relationship, then F(x) =e) mr DOr 
cases. So Z also has the marginal distribution F(x). The 


Beemer s a Consequence of the fact that m is defined to be 


ene Stationary vector associated with P. 


B. THE PRODUCT-MOMENT CORRELATION 


Theorem 2. 


The correlation coefficient between Y and Z becomes 





where 
=e B= A, 7 (l- a.) ds 
and 
_ _ 2 2 
M = (uy Uy) T1 1/0, ’ 
where 
sire 
F (x) 
u = f cl 
1 0 F(X) 


Wi 
UJ 





00 F(x) - F(x) 





nD ee eci(sc nl) aa 
x Oo 
Oo 
=O 
2 2 F(x) 2 
O = f Som Gl = hl 
l 0 F(x,) aE 
Deo c eon Ea) = F(xX)) 7 5 
29 l - F(x.) 2 
x Oo 
O 
o 2 = f[ x? aF(x) - EX] 
x 
0 
e 2 y : D 
= O71 ™) - O 5 T 5 7 (uy Uo) ™) To 
Proof 
; _ covlY,Z) 
YZ Oy on 


E(Y¥Z] - E[Y] E[Z] 


Ov On 


Now 
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apy 2 | 


E,(E(YZ|¥,Z « S]] 


= E[YZ|Y e X,,Z ¢ X, JPY ¢ X,,2 € Xj] 


alee 


TZ we X1,Z ¢ X,JP(YeX,,% € Xo] 


SED Za eee o oe X,IP(Y « X,,2 € Xj] 


oe 2 


Tey | Y ~e X5,Z ¢ X JPY ¢€ X,,2 € X,] 


= ELX, ] E[X,] Toy + E[X, ] E[X.,] my (L-a,) 


> E(X,] E[X,] m5 (1 - a5) + E(X,] E[X.] Ty Oo 


2 
Hy 4, Ty F By Uglloray) my 


2 
Fup Hgl(l-a5)t5 + uy” A575 


where 
Lee EIX,], Hy = E[X,] 
Further, 
miele = E_[E[Y|Y ¢ S]] 


Sethi Ye X,]P(¥ « X,] + Bay ee XIP(Y € X] 


= E[X, ] tT + E[X,] TT 


1 2 


TT = E[Z] 


Soll, liana uealic? 


by Theorem l. 


aD 





Also by Theorem l, 


2 


2 
ony ) 


= E[Y*] - (E[Y] 


If we put together these formulae into 


BIZ = EIY)E(Z] 


O — 
eZ Fy op 


we get 
(a, > me) Seen (Ch ie) 
meee) oe a 
and let 
ees Ca (1 - a.) 
Z ete 2 
M = (uy Uo) T1To/ Oy 


then 


C. GENERAL ALGORITHM 
We give here three algorithms for implementing the 


bivariate mixture-truncation method, which we call the Fxo 
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method, the UXO method, and the TXO method. All of these 
methods are exactly the same except in how the algorithm 
chooses Kos the truncation point, from the x, range [x x]. 
The first procedure, called the FXO, choose x, as a fixed 
moint of X 5 and x and uses the same Xo during the entire 
routine. The second procedure, called the UXO, chooses Xo 
mnicormly from [x 7X] and repeats this step in every call 
to the algorithm. The third procedure, called the TXO, is 
the same as the UXO procedure except in that it uses a tri- 
angular distribution instead of uniform. It is necessary to 
fix these choices of Xo becuase in general there is more than 


one X, which will give a bivariate pair (Y,Z) with the given 


marginal distribution and given correlation. The first 


procedure, FXO, is defective in terms of their discontinuity 


of distribution while the second and the third, UXO and TXO, 


are satisfactory in this respect. The choice of the midpoint 


of the interval Ix, 7x] OG XS in FXO is based on experience 


presented later for various marginal distributions. Note 
that the algorithm described here is inefficient in that it 


generates the truncated variables x) and X., by comparing ran- 


2 
dom variables X to x until one which is respectively greater 
than or less than x, is found. More efficient methods can 
be found in special cases such as the exponential, but the 
present algorithm requires only a generation of univariate 
random variables X without regard to the method used to do 


this. Of course initialization is required and this is 


Specific to each marginal distribution. 


an 





General Mixture-Truncation Method 
i Cimttlalization) 


4) For given marginal distribution F(x) and correla- 


2. Define truncation point Xo 
* FXO method 


| ns i 
1) oe x(x, +X) 


eo UXO method 


1) Generate a uniform [0,1] random variable V 


Eton COcCELIClent o find x, ranges [x 7X] 
1 


Bs ms A _ . 
11) Xe Xo (x, x, ) Vi 


* TXO method 


1) Generate two uniform [0,1] random variables 


ee 
11) a = X 5 + xy > Xo 
where 

ae oe 

Xn = 5 (x + x1) 
= = * 

x7 (x x, ) Vy 
= = * 

x5 (x, x1) V5 


3. Compute parameters value, Tyr Tor Ae 


4. Choose type for Y 
1) Generate a uniform [0,1] random variable U 
11) ieee. < 1 


1 9° cong 
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5. Y 1S an X, 
1) Generate a random variable X from F(x) 


mio tf X > a set Y « X and go to 6 


iii) Otherwise return to 5.1) 


6. Choose type for Z 
i) Set U+ ((U-1,)/(1-17))) 


ioe lt U < l-a GOeco. > 


oe 
mey. 24 1s an xX, 
1) Generate a random variable X from F(x) 


mye 6Lf X C> Xor set Z <« X and go to ll 


iii) Otherwise return to 7.1) 


Pewee2 LS On xy 
1) Generate a random variable X from F(x) 


moe 6 X < a set Z + X and go to ll 


111) Otherwise return to 8.i) 


fay 61S an xy 
1) Generate a random variable X from F(x) 


meee LL X < X,, set Y + X and go to 10 


fy) 6Otherwise return to 9.i) 


1Q. Choose type for Z 


i) Set U <« U/t, 
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mm) Shee < & go to 8 


I? 
iii) Otherwise go to 7 
ite Deliver (Y,2) and go to 4 for the FXO method, or go 


to 2 for the UXO and TXO method until a sufficient 


number of random vectors are obtained. 


me BIVARIATE DISTRIBUTION FUNCTIONS 


From Theorem 1, in Section III-A, we know that if Y 








is from Xi then 
BD (zy) See ey ce Lo) A 
F(z) 
4 Fix.) Met: 97 iS x 
O 
F(Z) - F(x) | 
yt Goo) TTF ete Z > X, 
maaeit Y is from oh then 
ez iy) = (l-a,) Fy (z) + a, Fx (2) 
i 2 
_ ECz,) : 
(l- a.) F(x.) if z< xX, 
Oo 
F(z) - F(x) 
(l- a.) orca aad Gh ize aa 
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By using these we can define the bivariate distribution 


ifunction as follows. 


Miy,Z) = PLY < y,@ < 2] 
Ng 
= fee ee ee) apy eu 
where 
me << ul) = Flu) 








a4 ra) If u< Xo, 
F(z) - F(x.) | 
Ger) en sey eas uU<X,, 
miz= 2|/Y =u] = ; 
(l= a,) ra equa 
F(z) - F(x) | 
(l- a.) + a5 T-FR) Te >Xoy 


So, if we put these together, integrating with respect to 


GP([Y < u], we get the final result: 
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a 
- RZ E CY) ates Y¥S%r ee (ITI-D-1) 
(1-a,) 
| fa, + 7 [F(z) -F(x,) J IF (y) Hee YSX%, 2>X, (III-D-2) 
(1-a,) 
PIY<y,Z<z] = ta, + 15 [F (y) -F(x,) ] F(z) rey >X1 ZS (ITI-D-3) 
mel 
a7, + (L-a,) [F(z)-F(x,)] a (III-D-4) 
Z 


a 
+{ (1-a,) + SUF (2) FC) J}(F(y)-F(x,)] if y>x,, 2> x, 


For example, the expression (III-D-4) is obtained as 


MY 
Ey ,Z) i Eee oe ly GE Cu) y > Xoe ZS x 


O 


ll 
ea 
rd 
N 

| A 
N 
K< 
il 
c 
A 


i xO] GE Gi) 


x 
a ecige 2 |v =u x] dF(u) 


il 
ro 
, 
+ 
ry 
N 

f 


F(x.) 1} A) 
a. 
' (l- a.) ts — = Eee ale ue (oy) = F(x,)] 


fe is easily seen from (III-D-2) that when z > o, 


mr < y-Z cee ys erom (I1II-=D-3) that as y + =, 


62 





P[Y < ~,2Z 


P[Y < ~,Z 


Peamecaromrrom (iti—-D-4) that as y > «, 


Z| 


cine | DZ meaneGmenat Zao, Ply < yi Zo< 17] = F(y). 


In particular from (III-D-4) we have that, as y > ~, 


my,z) > aj F (x) + (l-a,) [F(z) - F(x) ] + (l-a 


cc a, {F (2) ~ F(x) ] 


TT 


aj F(x) + (l-a 


2? Z 


F(z) + (l-a - (l-a 


2)" 


1771 


Seatee(OE(0 2p cay) 


Retz) 


2) 


where at the last step we used the facts that F(x) iy 


and m,(l-a)) = mo(l-a.). If F(x) is absolutely continuous 


with probability density function f(x) 


‘for the bivariate pair (Y,Z) is 
| 








a,(l-a,) | 
lea) de Gy;) fez) er Y 
L-a, 

= vee (2) Bet, Py, 

f{y,z) = 
L-a, 
Eye Cz) if y 

2 
— f(y) £(z) if y 

2 
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J A 


| A 


themrcne JOLNt p.daft: 


Crit =>) 


Ceri po) 


Cor) — 7) 


(iit >=) 





Note that there is a discontinuity in the density function 
as one crosses the boundaries of the four quadrants defined 
by the lines y = Xi Z =X. The density is the same in the 
first and third quadrants. The multipliers of f(y) f£(z)/r, 
are the same in all four quadrants iff there is independence. 


This occurs when (l-a,) = a5 SO that T) = and 


f(y,z) = f(y) f(z) for the whole range of y and z. 
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IV. BIVARIATE UNIFORM GENERATOR 


The uniform random variable is a continuous random 
variable with probability density function which is constant 
over the interval (a,b) and zero otherwise; the density 


function 


——— 3 Sf peat See Mi IS, 
boa —-_ = 
f(x) = 
GF otherwise ; 
mean = E[X] = a8, 
2 
variance = VAR[X] = Ps 3 


We note that if U has a uniform distribution over [0,1] 
then X = a + (b-a)U has a uniform distribution over [a,b]. 
SO we will only be concerned with algorithms for generating 
mattorm [0,1] distributions. 

Two algorithms for generating uniform bivariate pairs 
are given in Section II, Gaver's and Lawrance and Lewis's, 
and since they have relatively smooth distributed functions and 
are simple to generate, they are probably preferable to the 
uniform bivariate random variable generated by the mixture- 
truncation method unless Xo 1s taken to be random (of course 


mencorrelation for the Gaver bivariate uniform is difficult 
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to compute). However the development in this section of 

the uniform mixture-truncation method does help to fix ideas 
on the problem of initialization and determination of the 
possible range of correlations attainable in other cases. 

To use the mixture-truncation method for generating bivariate 
random vector (Y,Z) whose marginal distribution is uniform 
over [0,1] identically and has given correlation o, we should 
use uniform [0,1] distribution functions as F(x). Then Xi 
which is X constrained to be on the left side of the trunca- 
m2on point a becomes uniform over [0,x,] and X51 which is 
constrained to be on the right side of Xo becomes uniform 
over [x ,1]. 

for)|6 DETERMINATION OF PARAMETERS IN THE UNIFORM MIXTURE- 

| TRUNCATION METHOD 


Because xy 1S) un Gorm 10,x0] and x5 is umd form [x ,1], 


—— ay 





=o 
eS a pup ls 
ae 
Tanita ce = 
it 1 12 
Ll+x 
E(x, ] am aD = 2 
3 (1-x,)* 
VAR [X, | = G5 = 1D 
and 
l1- Ot 
pee oc’ ~ I-a,t1-a, ~ *o ’ eo 
1 2 
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= 1l- F(x.) — Loe = ier: 2 (IV-A-l1b) 


If we use these formulas in Theorem 2 of Section III, we get 


M = 3x (1-x_) 
Oo Oo 
a, - Xx 
2 ilk 
eo 1- x 
fe) 
op = 383 x Cay = x) (IV-A-2) 


Then from (IV-A-2), (IV-A-la) and (IV-A-1b) 


ae sta: Nes 
a eee. * Xo (IV-A-3a) 
Oo 
TT 
= _ 1)/1- 
to = 1 oe a4) 
Xx ox 
eo Oe a Ps (IV-A- 3b) 


1- x, 3(1- x) 


Furthermore, and a, are probabilities so they have to 


a. 2 


satisfy the following conditions: 


| A 
b= 


[A 
- 
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By solving these two inequality equations, we can get the 


feasible range of x for given p as follows. 


Lower bound of x 


* 
il 


Qi oO 
Vo V7 45 (1/3) 6 Icamiaom > 0 
es (IV-A-4) 
Yap fas if 7 omace() 
= Upper bound of x, 
1/2 + v1/4-(/3) 0 ihjaeoeee 2 
= (IV-A-5) 


Oyo ceo oO 


meter finding x, and x choose x, within the range (XX J 

by way of either the FXO, UXO, or TXO method. And by formulas 
(IV-A-la), (IV-A-1b), (IV-A-3a) and (IV-A-3b), we can get 

Tyr Tor Oy and Oo» Also we can expect some limitations in 


0 because of a From 


= 
it 


3 x, (1 - x) 
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We can compute that 


ee = by eee es = 1 x 


{A 
* 
il 


1/2 ee 


OD 
t 


7 2 


J A 
~*~ 
iI 

J A 
* 


In another way, note that if one chooses x, = it? 2 and 


Ao = 4, = 1, from (IV-A-2) we attain the highest allowable 


correlation p = 3/4, and when a, = 0, 9 = -3/4, the lowest 








attainable correlation. 
Memcurns out that for this choice of x, = 1/2 tin the FXO 


method), the bivariate uniform distribution is also about 





the smoothest, as can be seen ina later section (IV-D), 
| 


mma £rom the bivariate distribution. 


B. GENERATING PROCEDURE 

We developed here all of the three procedures, the FXO 
method, the UXO method, the TXO method for generating 
bivariate random vectors whose marginal distributions are 
uniform over [0,1] and correlation coefficient is p. In 
this algorithm we used direct generating procedures for X) 
and X, instead of comparing random variables X to XO uns 


One which is respectively greater than or less than x is 


found. 
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Uniform Mixture-Truncation Method 
(Initialization) 
1) For given -3/4 < 9 < 3/4, find x, and a 


lean eye = 173 p AE p > 0 


- = 

Y= "p73 if so) = 20 
- WeZ2etevi7e = 11/73 6 at » > 0 
u 


ian a3 if» < 0 


Define truncation point Xo 
* FXO method 


a) ee (x, +x )/2 


* UXO method 


1) Generate a uniform [0,1] random variable U, 


- _ o ‘5 
al) x, x, “= (x, x) Uy 


* TXO method 


1) Generate two uniform [0,1] random variables 


ya 2 
el. } Xo = x, ~ xy - x5 
where 
-_- a * 
xy (x, x, ) a 
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3. Compute Tyr Tor Oy and A. as 
Ta Gs F(x) = 36 
™T = US ™) 
Se 3 Xo i es 
a, = iL = fa - a4) 


4. Choose type for Y 


1) Generate a uniform [0,1] random variable U 


muerte U < nw. go to 9 


1 


fee y 1S on X, 


i) Generate a uniform [0,1] random variable Ww, 


oe _ ‘A 
li) Y< Xo + (1.0 x) Wi 


6. Choose type for Z 


i) Set U + ((U- m,)/(1 ~ ™))) 


ii) If U< 1- a,, go to 8 


ra 





ms. 2is an x, 


i) Generate a uniform [0,1] random variable W. 


A i) Z<+ Xo + (1.0 - x,) * Ww. 


fm) Go to ll 


fe 24 fis an x, 


i) Generate a uniform [0,1] random variable Ww. 


+ 4 * 
A 1 ) Z «+ XG W. 


1ii) Go to ll 


9. Y4i1s an xX) 


i) Generate a uniform [0,1] random variable Wi 
moe Y < x * Wi 


O 


10. Choose type for Z 


1) Set U < U/ 
mee if U < a, go to 8 
iii) Otherwise go to 7 
ll. Deliver (Y¥,Z) and go to 4, for the FXO method, or go 


co 2 for the UXO method and the TXO method until a 


sufficient number of random vectors are obtained. 


The programs are listed in the appendix and the scatter 
slot and bivariate histogram of resulting random vectors are 


shown in Section IV-D. 
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@e REGRESSION OF Z ON Y FOR GIVEN op 

Schmeiser (1979) in developing a bivariate gamma method 
fixed the free parameter by specifying the regression of Z 
and Y and vice-versa. We do this now for the uniform case 
with fixed Xo and uniformly distributed Xo: First in fixed 


x, Baoe we have, for y < Xo 


E(Z|¥=y,y< x] a, E(X,] + (L-o,) E[x,] 


i 
= 1/2 - =~ o 
6x, 


and for y > Xo 


E(Z|Y¥=y,y>x)] (1-a,) E[X,] + a, E(x] 


A a 
= CE TSO eee 
O 
These are the step function and clearly not linear in y. 
aT Xo has uniform distribution over [x , ex J, we have the 


following results. These results are dependent on Y values. 


me y < Xo then we have 
x 
u 
Pziy=y] = i EULA ce £(x,) qx. 
x 
ie 
a 
oe ee 1 
~ J (5 6X e) Xx -xX dx 
x, O u Q 


iS 





hale (a x HS x - 2 in 2B 


If xX, < y < x,, then we have 


ee 


by 
| E(Z|Y=y] = J E(zly=y, ee veenones (x Radice 
| x 
| L 
x 
u 
+ f E(Z|Y=y, eee eee cicioc 
ye 
¥ p 
eo cue) *\%,) 2 
x, 
x 
oO 4 ’ 
+ f Cree) £ (x) dx 
y O 
1 ip 0 ey. u 
= = (= x = x -= in ==) 
x7 %, 2 2 ye ok x 


im: y > xs then we 


x 
u 
E(Z|Y=y] = rene age ey cel) (acm) doe 
x 
L 
x 
~ J Emcee 5) 1° %o 
MD 
ba Xx 
= 1 1, iy _2 u 
me wnt > 2 Gg -! eae 
u R x 





‘From the results we can see that if we use uniform distri- 
ition for x we can get smoother regression functions if 
although we still have step functions for the 


ee = *y 


y < X, and y > xX, cases. 


D. SIMULATION RESULTS 

We will show here the scatter plots and the bivariate 
histograms of the mixture-truncation bivariate uniform (0,1) 
random vectors with correlation op = 0.3 and p = -0.3 by the 
FXO, UXO and TXO methods in Figures (IV-a), (IV-b) and 
(IV-c), respectively. 

In addition, we will show the scatter plot and bivariate 
histogram of sample size 2000 from Gaver's transformation and 
Lawrance-Lewis method in Figures (IV-e) and (IV-d), 
respectively. The Gaver method does not have known correla- 
tion but the value of p to give op in the exponential case 
1s used. 

As we mentioned earlier, the FXO method has some discon- 
tinuity at the truncation point. But the UXO and TXO methods 
generate relatively smooth continuous distributions. In this 
respect, we say that the FXO is defective and the UXO and 
the TXO are relatively satisfactory. 

The computed correlation from subroutine BIVHST (Bivariate 
histogram) is a little different from given correlation. But 
this can be assumed as a sampling error. To check this error, 
we simulated 10 times with given correlation p = 0.1. From 


this simulation we get the following results: 
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= mean of computed correlation = 0.105 


0 | 


VAR[p] = variance of computed correlation = 0.0004 


cofo] = standard deviation of mean = 0.0064 


uo 








Figure IV-al. Scatter plot for sample of size 2000 from 
mixture-truncation bivariate uniform with 
0 = 0.3. Here x_ is fixed at the midpoint 
between the lower and upper bounds. 
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Scatter plot for sample of size 2000 from 
mMixture-truncation bivariate uniform with 
0 = 0.3. Here x_ is fixed at the midpoint 
between the lower¥ and the upper bounds. 
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Figure IV-bl. Scatter plot for sample of size 2000 from 
mixture-truncation bivariate uniform with 
op = 0.3. Here x. is taken to be uniformly 
distributed betwéen the lower and the 
upper bounds. 
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Figure IV-b3. 


Scatter plot for sample of size 2000 from 
mixture-truncation bivariate uniform with 

0 = -0.3. Here x_ is taken to be uniformly 
distributed betweén the lower and the upper 
bounds 
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Segure IV-c3. Scatter plot for sample of size 2000 from 
mixture-truncation bivariate uniform with 
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V. BIVARIATE EXPONENTIAL GENERATOR 


Another probability distribution of major interest in 
simulations is the exponential. The cumulative distribution 
function and probability density function for the exponential 


are, respectively, 


F(x) 


il 
t= 

{ 

(D 

*% 

|v 
© 


and 


£ (x) Tok 


H 
> 
D 
~ 

[v 
© 


The expected value of the exponential distribution is 


} 
H 


Etx}) >= I1/r 


and the variance is 


The problem of generating exponential deviates reduces to 
@ee,OL generating “unit" exponentials, i.e., those with 

A= 1.0, and then multiplying the result by whichever 1/) 
is necessary to give the desired distribution. That is, if 


the random variable E has the exponential (i = 1) distribution 


then, 
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nas the exponential (A = n) distribution. In this section, 
we conSidered only unit exponentials as marginal distributions 
for bivariate pairs. 
eee DETERMINATION OF PARAMETERS IN THE EXPONENTIAL MIXTURE- 
TRUNCATION METHOD 
Because xX) 1S an exponential (A = 1) truncated to the 
Hert of Ns and Xo is also exponential (i = 1) truncated to 


Sie right of Xo, we have 











x 
O 
= a F(x) 
BIX i = < * q F(x_) 
TT 
— ae ee ’ 
T, “0 
Xo 
2 2 F (x) 
VAR[X,] = a noe ex id u 2 
i: dL 0 F(x,) l 
T 
. he 2 
= 1] 5 Xo ; 
oat 
co F(x) - F(x) 
cere eo, eri 
x O 
O 
= Jl+ xX; 
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VAR [X, 


mee from definition, 


oat 


lI 


7 2 2 
ee seri YD 
Xx O 
Oo 
= ie 
1 - a, 
l-a,t+tl-a, 
F(x) 
Le P (V-A-1la) 
L- a, 
L-a,+1l-a, 
1- F(xX,) 
e : a = AS)) 


If we use these formulas in Theorem 2 in Section III, we get 


and 


~ 
se x A (V-A-2) 
TT © 

nf 

a 
—(a 1T. ) (V-A-3) 
Ts ik iE 


2) 





= x (V-A-4) 


For given correlation coefficient o, we can compute %, as 


@eeeuniction of Xo from the formula (V-A-4) as 





OT 
i 
e% = ate ah 
1 me 1 
O 
= (l+—y)(i-e °*) 
x 
O 


T 


and we know that a. = 1-==(1- a) from the formulas (V-A-1a) 


2 
and (V-A-lb). From this, 


TT 
1 - —(1-a 
2 


a5 1? 


These 04 and a, are probabilities, so they have to be greater 


than or equal to zero and less than or equal to one; 
eoieese -)(l + —~) « 1 (V-A-5a) 
Memee cater 6 °)° 8. <« 1 (V-A-5b) 

From these two inequality equations, (V-A-5a) and (V-A-5b), 


= Can find the x, ranges for given correlation coefficient 


Pp. To solve these equations, we can divide into two cases, 
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one for positive correlation case and the other for negative 
correlation case. If correlation coefficient p is positive, 
then both equations are always positive. Thus we only need 


to find the x. ranges which makes (V-A-5a) and (V-A-5b) less 


O 


than 1. From the equation (V-A-5a), for the a, case, 
rio 
a, = (l-e )(1L + =p) See 
x 
O 
becomes 
Xx 
O 2 
ole IL) eg x, 
and let 
x 
_ 2 _ Oe 
oa) os —o ’ MS = pe 1) 


Because of the first derivatives of Y1 and Y5 are always 
positive, we know that these two functions are monotone 
increasing functions. Thus we can find x, ranges which 
Satisfy Wy 2 Y> by the Newton Raphson method. When using 
the Newton Raphson method, let y = ie: = QO and find an 
approximate solution, by approximating exponential series, 


which we can use as a starting point. That is, 


2 3 
2 ac 6 
Meee DCC CME XO Porta 7 
mee, 2 _ Bo) = 
= ¢ XS (1 >) Xo +p = Q 


oa 





Then 


eS ea C= = ae 


O WW Sep 


Starting with this approximate value in the Newton Raphson 


method, we can find X, range, say [x, rX }, which satisfies 
i i 


O<a, sl. And for the a, case, 


LA 


becomes 


This result is exactly the same as the a, case, that is, at 


ul 


the same range Oy and a, satisfy constraints A, 1 and 


2 
a, < Ll. ThusS we can use x and x as the lower bound of 
oo'*," and the upper bound of XorX ys If correlation coeffi- 


cient p is negative, then the equations (V-A-5a) and (V-A-5b) 
are always less than 1. Therefore we need to consider only 


one constraint which makes a, > 0, a5 > O. From the a, 


equation (V-A-5a), solve the inequality equation 


—Xx 
cme ee. 
x 
O 


Oe ay 
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) 
=e 
O 


| Mince = ¢ is always positive, we see that to satisfy 


the inequality l + a should be positive, l.e., 


x 
O 


~—; > -] 
x 
© 


equivalently, we have 


In the A, Case, from equation (V-A-5b), 

ae ve ~aS Ze 
is A, = e +e “(l-e ) a} 
xX 
O 
or, equivalently, we have 
x 
x - . mee 


As in the positive correlation case, we can find a starting 
point by approximation to solve this equation by Newton 


Raphson. The result comes out as 


ryan o)/ (= 5) 


Ss 


with this starting point we find another bound of Xo which 


Satisfies 0 < a This becomes the upper bound of Xr Xe 


Ze u 
and from the a, case, we have a constraint eee which 
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becomes the lower bound of Xo ee eer ee Vv-o. The lowest 
and highest correlations available for bivariate exponential 
pairs in mixture-truncation method are approximately -0.480 and 


0.647 respectively. By comparison note that the most negative 


correlation available for bivariate exponential pairs with 
2 


identical fixed marginals is ers [Moran (1967)]. Gaver's 


(1972) negatively correlated pair has correlations in the 
range [-0.5,0]. The table (V-1) shows the lower and upper 
pound of x in the mixture-truncation method with identical 


marginal exponential and given correlation. 


Table V-l: The lower bound and upper bound of x 
in the mixture-truncation method witR 
identical exponential marginal distributions 


x, range for each correlation 





B. GENERATING PROCEDURE 

The generating procedure of bivariate exponential random 
vector is almost the same as in the uniform case. As in 
the uniform case, we develop three procedures which are the 


FXO method, the UXO method and the TXO method. As mentioned 
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in Section III, all of these methods are exactly the same 
except in how x is chosen from the x, range [x,,x,1- The 
FXO method chooses x, as a fixed point between x, and x 
This procedure makes some discontinuity at the truncation 
point. In this respect, this method is defective. 
The UXO and TXO methods are the same except we use 

different distributions for a In the UXO method we use 
the uniform distribution and the ERianogUlar -GiStrabution for 


the TXO method. These methods have more smooth distribution 


than the FXO method. 


PePonentlal Mixture-Truncation Method 


i. (Initialization) 


ie Bor given -0.48 < 9 < 0.64, find x, and x 


2. Define truncation point XS 
* FXO method 


i) ns = (x + x ) 


* UXO method 


i) Generate a uniform [0,1] random variable Uy 


ee, _ _ e 
ea: ) Xo x, + (x x, ) U, 


* TXO method 


1) Generate two uniform [0,1] random variables 


Var Vv. 
11) Xo = x, + xy + X5 
where 


1Bjat 





Xx = (x, + x) /2 


= (x, =e V 


3. Compute parameter values 


™) = F(x) = l-e 
T > = l- ™) 
a) = m, (1 te a 
Oo 
TT 
OL = l- a aay) 
2 To i 


4. Choose type for Y 


1) Generate a uniform [0,1] random variable U 


mm If U < go to 9 


1 
Be Y iS an X, 


1) Generate an exponential random variable EL 


ii 6CL E E) > Xr set Y + E) anamce CO 6 


1ii) Otherwise, return to 5.i) 


6. Choose type for Z 


i) Set U <« ((U-m,)/(1-7))) 


me If U < 1 =a Ont Oma 


ps 


f 2is an X. 


i) Generate an exponential random variable E. 


may 6|6CL EL E. > Xo set Z + E. andegO tO ji 


iii) Otherwise return to 7.1) 


8. Zi1s an xy 


i) Generate an exponential random variable E, 


ii) If E~. < x, set Z+«eE 


9 £ % andsdGuco Li 


2 


mer Otherwise return to 8.1) 


i DG 
we Y 1S Aan Ay 


1) Generate an exponential random variable Ey 


, s@Ct Y « E 


a anGad@ te 1.0 


11) ete Ey =x l 


111) Otherwise return to 9.1) 


10. Choose type for Z 


i) Set U < U/T, 


marr U < a,, go to 8 


Le 


111) Otherwise go to 7 


ll. Deliver (Y,Z) and go to 4 for the FXO method, or go 
to 2 for the UXO and TXO methods until a sufficient 


number of random vectors are obtained. 


For the exponential case it is possible to give a more effi- 


Cient algorithm in which Ay and X, are generated exactly. 


7 
The algorithm is as follows. 
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Efficient Exponential Mixture-Truncation Method 


o (Initialization) 


i) For given -0.48 < p < 0.64, find X and x 


2. Define the truncation point Xo 


* FXO method 
ae 
Xo = ZX, ata x.) 


* UXO method 


1) Generate a uniform [0,1] random variable U, 


im) xX = x + (x - x ) * U 


Oo Q Q iL 
* TXO method 


1) Generate two uniform [0,1] random variables 


ane} 
11) Xo = x, + x5 == x5 
where 
x, = (x, + x) 72 
z » * 
x4 (x, x, ) Vy 
cs zs * 
Xo (x, x.) Vv. 
3. Compute parameter values 
™) = F(X) = l-e 
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2 nt 
0 
Xo 
TT 
a = l- ae: - a.) 
2 T 5 1 


Choose type for Y 


i) Generate a uniform [0,1] random variable U 


ii) If U< 7), go to 9 


Y ais an x5 


1) Generate an exponential random variable Ej 


11) Set Y « x +E 
Oo a 


Choose type for Z 


ie cet U « ((U-m,)/(l-1 )) 


i 


ieee tf U < 1 - o Goviro 


et 
21s an x, 
i) Generate an exponential random variable E. 


1i) Set Z < x, + Es ancmgorc@ i 


Zi1S an x) 


1) Generate a uniform [0,1] random variable Ww. 


11) Set Z@<« = 1ln(1.0 - W,*t,) ania GO tO. LL 


Y 1s an xX) 


1) Generate a uniform [0,1] random variable Wi 


05 





Pome oete Ze — intl. 0 = W. Xn 


10. Choose type for Z 


1) Set U < U/ny 


mae if U < «a Tome oo 


ee 


iii) Otherwise, go to 7 


ll. Deliver (Y,Z) and go to 4 for the FXO method, or go 
to 2 for the UXO and TXO methods until a sufficient 


number of random vectors are obtained. 


Note that to compute x, and x, iste pulesot both algorithms 
we use subroutine BOUND which is listed in the appendix and 
the program with effective algorithm also. The scatter plot 


and bivariate histogram of resulting random vectors are shown 


in Section V-D. 


m REGRESSION OF Z ON Y FOR GIVEN o 

Schmeiser (1979) has used the regression of Z on Y=y 
to fix the parameters in his bivariate gamma distribution. 
‘Consequently we investigate this for the mixture-truncation 
method case. The regression is different depending on 
meether Y < x, Or Y > Xo- We consider two cases here, one 


for fixed x and the other for x Dapeng sn LfOrm- adiLStribution. 


For fixed Xo we have 


E(Z|Y¥=y,y<x J = a, Elx,] + (1-«4,)E[x,] 


TT 
Z 
x (1 + =e) 


= lt Xx, 7 4 
ak 


1 


106 





Substituting the value for a 


a 
_ O 
a, = mm (1 + 2 
O 
then we have 
_ = = p 
E(Z2|Y¥=y,ysx,] = 1+ x x, (1 +85) 
x 
O 
= Ll saat 


mnad if y > xor then 


E(Z|¥=y,y>x)] = (1-o,)E[x,] + a,E[x,] 
TT x 
= 1-4 x, + 2a, 
1 mt 
Substituting the value for tos 
TT 
= ee 
A, = i T° a5) 
Tl 
= . Bio OS: 
= le m (1 _ =) 
x 
O 
then we have 
| a 
men Y—yey>x ) = J] + — —- 
O T> Xo 


HO 





Thus the regression is constant over (0,x,) and changes for 
=> X.- This is not surprising in light of the joint dis- 

mribution given in Section IV-D. For uniformly distributed 
a the computation is different for different ranges of Y. 


me y < x,, then we have 





ee 
x 
u 
migiyY=y) = era av X=xX, Y <x] ceo dx, 
Xx 
g 
x 
: 0 6 
— s (1 = menor S ax, 
Xo Xo 
-X -X a rc 
g 
= - e ae eae f = xy AX, 
Xp a 
-x -X -X —X 
= e Se Ghee yt Aaa ea ae ge i 
x xX 
u g 
4 
+ ln —- x +X) 
x i c 
g 
fc 
eo 


then we have 
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ye 
miz|Y=y) = ME =, X=x,, Y>2x ]£(x,) dx, 


*y, 
x 
u 
+ f E[Z|Y=y, X=Xo, YS xX J£(x,) dx, 
Mf 
xX 
Es =v 7 os te Ley: 
= 2e + 2py eS. e y~ ) 
*u 
+ Ue 
' y 
ma yy > x then we have 
x 
u 
mie1Y=y) = s PigyY = ye X= xX, y> x J£(x)) dx. 
Oe 
Q 
a T —X 
Geri = 4) a. ° ax 
ite oO 
x Zane 
i 
“xX, iG 
aCe E oeitc Gy.) 


By making a uniformly distributed over the available range 
Gr x, for given correlation, we can get smoother behavior 


for the regression function. 


D. SIMULATION RESULTS 
We will show here the scatter plots and bivariate 


historgram of the mixture-truncation bivariate exponential 
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memcom Vectors with correlation p = 0.3 and p = -0.3 by 

the FXO, UXO and TXO methods in Figures (V-a), (V-b) and (IV-c), 
respectively. Also we show the results from Gaver's bivariate 
exponential and Marshall and Olkin's in Figures (V-d) and 

(V=e), respectively. 

In the mixture-truncation bivariate exponential distribu- 
tion, the generated random vectors by the FXO method also 
has discontinuity at truncation point but the other cases 
have relatively smooth distributions. 

The computed correlation printed in the left side of the 
bivariate histogram is a little different from the given 
correlation. But this difference can be assumed as a 
sampling error. 

To check this error we simulated 10 times with given 


correlation p= 0.1. From these simulations we get: 


p = mean of computed correlation = 0.092 
VAR[p] = variance of computed correlation = 0.0008 
pce = standard deviation of mean = 0.009. 


iO 
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Scatter plot for sample of size 2000 from 


Pagure V-al. 
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Here x 
between the lower and the upper bounds. 
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Figure V-a4. Bivariate histogram for sample of size 
2000 from mixture-truncation bivariate 
exponential with p = -0.3. Here x is 
fixed at the midpoint between the Yower 
and the upper bounds, 
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Scatter plot for sample of size 2000 from 
mixture-truncation bivariate exponential 
with 9 = 0.3. Here x_ is taken to be 
TAAAlein@ ta Tey’ distributed°between the lower 
and the upper bounds. 
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Figure V-b2. 
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Figure V-b3. 





Scatter plot for sample of size 2000 from 
mixture-truncation bivariate exponential 
with 9 = -0.3. Here x_ is taken to be 
uniformly distributed Between the lower 
and the upper bounds. 
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Scatter plot for sample of size 2000 from 
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with 9 = 0.3. Here x_ has triangular 
distribution between fhe lower and the 
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VI. BIVARIATE GAMMA GENERATOR 


The gamma distribution with shape parameter r and 


scale parameter A has the density function 


f(x) = i veESO= dl") 


meoecumulative distribution function 


=X 


F(y) ——- xX e ax. (VI-0-2a) 


It 
—_~ 
> 
Ky 
i 
~ 


In particular, if the shape parameter r is a positive integer 


then 


te, AY J 
ae eee ie oat (VI-0-2b) 


2° 


where [(r) is Euler's gamma function 


meee ne random variable X has Gamma (r,i) distribution then 
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when r = l, X has the exponential distribution while X, 
Suitably scaled, has an asymptotically normal distribution 
aS r > ~-, We note that if X has a Gamma(r,l) distribution 
then = has a Gamma(r,A) distribution. So we may set = 1 
as far as the generating algorithm is concerned. The output 
from the generator may then be appropriately scaled. Two 
algorithms for generating bivariate gamma random vectors are 
given in Section II. Since they have complicated computa- 
tions to define the parameter values and use inverse trans- 
formation functions, the mixture-truncation method is 
probably preferable to those methods. To generate bivariate 
random vectors with identical Gamma (r,1l) marginal distribu- 
tions and having given correlation p by mixture-truncation 
method, we only need univariate gamme (r,l) random variate 
generator. Under the assumption that the univariate gamma 
generator is available we developed this procedure. In fact 
we used Lewis and Robinson's gamma generator as a univariate 
gamma generator. 
A. DETERMINATION OF PARAMETERS IN THE GAMMA 

MIXTURE-TRUNCATION METHOD 

If the marginal distribution is gamma with positive 


integer shape parameter r and scale parameter ’ = 1, then 
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we know that F(x) becomes cumulative distribution function 


of gamma (r,l). Then the xy 1s gamma (r,l) truncated to the 


Bett of x, and Xo 1s also gamma (r,l) truncated to the right 


of Xo: Thus we have 





E[X, ] = 


| 
— 
~ 
OF 


a 
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E[X,] 


Let 


R2 


S 


—xX 


—~X 


O 


l|lroh 


r+1 


} 


1=0 


Moric 


(mai). 6 


(r+l1) ! 


(r+1-1) ! 


3.0 


Ea 


we! r-1 
(r-1)! “o 


x) - F(x.) 


1 - F(xX)) 


r-1 


O 


r+1 (r+1)! 
Oo ——————— 


ixo (r+1l—-1)! 


r+l-1 


O 


r+l-i 





and substitute, then we have 


2 ar es 
owe = ey 10 (x) : 
VAR[X.] = oa 2 — Artl)! - B_ (rr! - A)? 
1 1 mf (x) T r(x) 2 
= e A 
ee 
2 
Z B A 
VAR [X ] = O = ee a ; 
2 2 woP (xr) nT (x) * 


If we use these formulae in Theorem 2 in Section III then 


we have 
(ror! = a)? 
M = i... 
T 75 (T(r) ) ie 
and 
a. —- 7 
= a l 
8 ae “+ w 
2 
Thus 
(a4 - T,) (ror! - aA)? 
a Cr (VI-A-1) 
x ree 
T™] 5 
Or given o, x, and r, 
Oe ire rel) | Ty) T° 
oa = AF es (VI-A-2a) 
(7. le A) 


eS 





and 


oO = C7 R215) 


The value of (VI-A-2a) and (VI-A-2b) are determined with 
given o, x andes tnus £Or given Marginal distribution and 
correlation coefficient p we can find xX ranges as before. 


For simplicity, we will consider them with marginal distri- 


mition gamma (2,1). From the formula of (IV-0-2b), 
irs iG 
™ = F(x,) = l-e - X¥5 e 
-X ~X, 

To = 1 - ™ = = fe Xo e 

and 
“Xo 2 
A = e (x + 2x + 2) 
O O 


Then (VI-A-2a) and (VI-A-2b) become 


Z 
2 50 ™) TT 
eT = : =e = Ty) (VI-A-3a) 
x e 
O 
Zap r,° T 
O5 = T + ; Se . (VI-A-3b) 
x = 
O 


These Os and a, are probabilities, so they have to be greater 


than or equal to zero and less than or equal to one; 


P32 





_ 1 
OS 9% 5 sie FT St 
Soa 
O 

Zo n° T 5 
QO < Gy = Ts + ; Ene ae 

en 

O 


From these two inequality equations we can find the Xo 
mes LOL given correlation coefficient p. If p > 0, then 


ay and a, are always positive. Thus we only need to find the 


xX, ranges which makes (VI-A-3a) and (VI-A-3b) less than one. 


mrom equation (VI-A-3a), the a, case, and from equation 
(VI-A-3b), the a, case, we have exactly the same inequality 


equation; 


x 


Oo 
2pe (1 +x) < Xo + 20 (l+x)) (VI-A-4) 


Meee, «Chen a. and a. are always less than one. Thus we 


iL 2 
only need to find the x, ranges which makes (VI-A-3a) and 


(VI-A-3b) greater than zero. From the Oy equation, we have 


Vion 4 x) aX (VI-A-5a) 


where 


| o | 


and from the 5 equation, we have 


13:3 


Ti ae x + v2p*(1 + x.) (VI-A-5b) 


where 


Note that if we let ia for the left side equation of inequality 
Sign and Y>5 for the right side equation of inequality equa- 
tions (VI-A-4), (VI-A-5a) and (VI-A-5b), then we know that all 
Yy and Y5. are monotone iiereaisaind £unectiens. Thus as in the 
exponential case we can find the range of x by the Newton 
Raphson method. For the positive cirrelation case, equation 
(VI-A-4) will give the X, and x and for the negative corre- 
lation case, equation (VI-A-5a) gives Xo and equation 
(VI-A-5b) gives xX By using subroutine GBOUND which is 
listed in the appendix, we can find Xo and x: The tables 
(VI-a) and (VI-b) show the results from subroutine GBOUND. 
Unfortunately we have limitations on the possible range of 
correlations; we can generate bivariate random vectors with 
Meese tations in approximately the range -.5 < o < .6 for the 
Gamma (2,1) marginal distribution. In the same way as the 
above, the values x and x, can be computed with increasing 
complexity for integer values of r greater than 2. In 
principle it is also possible, using numerical integration, 
to compute allowable values for non-integer values of r. The 
computation is not as complicated as those required for 
schmeiser's bivariate gamma. Moreover the mixture-truncation 
method does not require that the user be able to compute 


the inverse gamma distribution function. 
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Table VI-a: x, ranges with Gamma (2,1) 
marginal distribution 
and given 0 


u 
3.41 
2S 
2.34 
Za 
en 





Table VI-b: x ranges with Gamma (3,1) 
marginal CrSeributivoen 


and given o 





oD 





B. GENERATING PROCEDURE 

The generating procedure of bivariate gamma random vec- 
tors is almost the same as the uniform and exponential cases. 
We also developed here three methods which are the FXO 
method, the UXO method and the TXO method. As in the uniform 
and exponential cases, these three methods are the same 
except in the choice of X, from the x, range [x,,x]- The 
FXO method chooses x, as a fixed point from the x, range and 
uses it during the entire routine while the UXO and TXO 
methods choose another Xo in every routine. From experience, 
the midpoint of the x range gives the best result of the 
FXO method. In the UXO method and the TXO method, we assume 
that x, has the uniform distribution and triangular distri- 
bution, respectively, over [x ,/x]. This assumption removes 
some discontinuity which occurs at the truncation point in 
the FXO method. 

Gamma Mixture-Truncation Method 
(for integer shape parameters) 
fee (initialization) 
1) For given allowable correlation coefficient 0, 


find x. and x 
g u 


2. Define truncation point x 
* FXO method 
. eee 
on) | xX = 5 (x, + x,,) 


O 


3 © 





* UXO method 


i) 
zd ) 
mee XO 


a) 


ah) 


Generate a uniform [0,1] random variable U, 


= = * 
x, x, 7 (x x) Us 


method 


Generate two uniform [0,1] random variables, 


Vie V5 
x, = x, + x1 - X5 
where 
x, = 1S tr x) 72.0 ; 
= Pah * 
x) (x x) Vi ; 
= = * 
x5 (x x) V5 


Compute parameter values 


where 


. e re 
TT a cea | —— , 
1 O 5=0 ane 
TH = i= T 
; _ One leas 1)! TT) 7 a 
1 (nox! - a) 1’ 
TT 
a = l- = Bian) 
T 5 1 


i / 





(r-1)! "oO 


4. Choose type for Y 


i) Generate a uniform [0,1] random variable U 


mete U =< 1), go to 9 


ee Y is an x, 
i) Generate a gamma(r,l) random variable G) 


ioe }6|6LL ft G) > Xs set Y + G ane egOuto 6 


iii) Otherwise return to 5.1) 


6. Choose type for Z 


me set U < ((U-7,)/(L-7,)) 


meet U < l-a,, go to 8 


a 


tee 2 is an X. 


1) Generate a gamma (r,1l) random variable G, 


a) CULT G. > Xo set Z < G. stare, je) weler JUL 


111) Otherwise return to 7.1) 


See 1S an Xx 


1) Generate a gamma (r,l) random variable G. 


moe Lt G, S Xos set Z < G. anaagourte. £E 


111) Otherwise return to 8.i) 


IL ie) 


ee Y iS an Xy 


i) Generate a gamma (r,1l) random variable G, 


71) ee Co 7 set Y <« G 


l emiGme cote pn” 


dL 


1ii) Otherwise return to 9.1) 


10. Choose type for Z 
i) Set U < U/r4 


fie, If U < « go to 8 


ee 


iii) Otherwise go to 7 


fie Deliver (Y,Z) and go to 4 for the FXO method, or go 
to 2 for the UXO and TXO methods until a sufficient 


number of random vectors are obtained. 


For the non-integer shape parameter case, step 1 and step 

3 need more complicated computation procedures but here we 

are only concerned for integer shape parameter cases. In 
step 1 we used the subroutine GBOUND which is listed in the 

appendix. The scatter plots and bivariate histograms of the 

resulting random vectors from this algorithm are shown in 


the next section. 


Se. oOoIMULATION RESULTS 

The scatter plot and bivariate histograms for random 
vectors of the mixture-truncation bivariate gamma variables 
with marginal gamma (2,1) distribution and given correlation 


(0.3 and -0.3) are given here. As in the uniform and 
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exponential cases, the UXO, TXO and FXO methods are used 

and are shown in Figures (VIF-a), (VI-b) and (VI-c), respectively. 
The FXO method has some discontinuity at the truncation 

point, but the UXO and TXO methods appear to give a relatively 

smooth continuous distribution. The computed correlation in 

subroutine BIVHST (Bivariate histogram) is a little differ- 

ent from the given correlation. But we can assume this 


difference is a result of sampling error. 
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Vite CONCEUS FON 


The mixture-truncation method is a general method which 
can generate bivariate random vectors having any theoretical 
Marginal distribution and allowable correlation. The 
generating procedure is very simple and doesn't need much 
computation for defining parameter values. In this respect, 
the mixture-truncation method is a very attractive method for 
generating bivariate random vectors. A price is paid for 
this simplicity and generality in that the Frechet bounds of 
correlation for the bivariate distributions specified by the 
Marginal distribution given by Moran (1967) are not always 
attained. Also there is some discontinuity in the bivariate 
distribution. However this discontinuity can be decreased 
by giving some distribution to the truncation point over 
its range for given p. Thus the mixture-truncation method 
is very attractive for simulation studies involving only 
partly specified dependency structures. The mixture- 
truncation method may be extended to generate bivariate 
random vectors having negative values. Another extension 
may be made to use grade correlation or rank correlation 
which are invariant under transformation instead of using 


the product moment correlation. 
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